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Abstract 

We address analytically and numerically the problem of crack path prediction in the 
model system of a crack propagating under thermal loading. We show that one can explain 
the instability from a straight to a wavy crack propagation by using only the principle of 
local symmetry and the Griffith criterion. We then argue that the calculations of the stress 
intensity factors can be combined with the standard crack propagation criteria to obtain 
the evolution equation for the crack tip within any loading configuration. The theoretical 
results of the thermal crack problem agree with the numerical simulations we performed 
using a phase field model. Moreover, it turns out that the phase-field model allows to clarify 
the nature of the transition between straight and oscillatory cracks which is shown to be 
supercritical. 



1 Introduction 



Crack path prediction is one of the main challenges in the field of fracture mechanics. The reason 
is simple - a satisfactory equation of motion of a cra ck tip is associated with a fundamental under- 
standing of material separation mechanisms (Freun d. 1990; Broberg . 19991 : Fineberg and Harder . 
19991 : lAdda-Bedia et al.l . Il999l : iLeblondl . 120031 ) . From a more general point of view and alongside 
its importance in many applications, the determination of crack propagation laws is necessary 
to describe the fracture phenomenon as a pattern formation process induced by mechanical 
stresses. Within the framework of Linear Elastic Fracture Hechanics (LEFH), the propagation 
of a crack is maiii l y governed by t he singula,r beh avior of the stress field in the vicinity of its 
tip dFreundl . Il990l : lBrobeT3 . Il999l : ILeblondl . l2003l h For a two-dimensional quasi-static crack, 
which is the main purpose of the present work, this behavior is given by 



(1) 



where Tilj{(f>) and are universal functions describing the angular variation of the stress 

field, and Ki and Ku are the stress inte nsity factors (SIFs). The evolution of the crack tip is 
gover ned by the Griffith energy criterion (jCriffithl . Il92d : iFreundl . Il990l : iBrobersJ . Il999l : ILeblondl . 
2003), which states that the intensity of the loading necessary to induce propagation is given by 
G = T, where G is the energy release rate and F is the fracture energy of the material, that is 
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the energy needed to create new free surfaces. This criterion can be rewri t ten u sing the stress 
intensity factor, in which case it is referred to as the Irwin criterion ( Irwin . 19571 ) 



Ki = Kic = y%ir 



(2) 



where Kjc is the toughness of the material and a is the shear modulus. The Griffith criterion was 
originally formulated for quasi-sta tic crack propagation (Griffitl^, M), and later generahzed 
to rapidly moving cracks ( Freund . 199d ). While this criterion is very useful in predicting crack 
initiation, it cannot predict the direction of the crack tip, and therefore in most cases it is not 
sufficient to determine the actual path of the crack. In order to achieve this, several suggestions 
have been made. Among them, the Principle of Local Symmetry (PLS) states that the crack 
advances in such a way that in-plane shear stress vanishes in the vicinity of the crack tip, or 
explicitly 

Kii = 0, (3) 

This rule was proposed for in-plane quasi-static cracks (iGol'dstein and Sakanik". ' 19741 : iLeblondl . 

and then generalized to rapidly moving cracks (lAdda-Bedia et al. . J_999|). Note that 
from a historical point of view, the first criterion for crack path selection, based on symmetry 
arguments similar to the PLS, was first formulated for antiplane (mode III) crack propagation 
(jBarenblatt and Cherepanovi . Il96ll ) . Recently, the dynamics of a rough crack in t wo dimensions 



has be en successfully described by an equation of motion derived using the PLS (jKatzav et al 
Another suggestion, based on symmetry a rguments, and which recove rs the principle of 



local symmetry in a certain limit, was proposed by Hodgdon and Sethna ( 19931 ). who formulated 
an equation of motion of the crack tip. However, this formulation introduces additional length 
scales, which are not known a priori. A di fferent approach, known as the maximum energy 
release rate criterion ( Erdogan and Sih . 19631 ) states that the crack advances in a direction that 
maximizes its energy release rate. Interestingly, even though this approach is quite different 
from the PLS, it produces very similar results t o those obtained using the PLS to such an extent 
tha t they were even conjecture d to coincide (jBilbv and Gardewl. Il975l). Howe ver, as shown 
in ( Amestoy and Leblondl . 19921 ) for kinked cracks and in dKatzav et Zl . [20073) for branched 
cracks, the results are not exactly the same and a clear distinction can be made bet ween the 
two. Actu ally, it has been shown that the PLS is the only self-consistent criterion (jLeblondl . 
1989 . 2OO3I ). A theoretical challenge is to explain the current rich crack path phenomenology, 
including the various known instabilities, using a pure LEFM approach combined with the 
Griffith criterion and the PLS. This has been recently argued for roug hening instabilities o f 



cracks in dKatzav et al.1 . l2007al ') . as well as for the branching instability in dKatzav et al.1 . l2007bl ) 
and here we wish to contribute further to this effort. 

Trying to address the problem of crack path prediction from the experimental side, the ther- 
mal crack problem has attracted a lot of attention. The propagation of cracks induced by ther- 



mal gradients has been widely studied (iRonsin et a .1.11995 : Ronsin . 19961: Yuse and Sano j_ 1997 
Ronsin and Perrin . 1998 : Yang and Ravi-Chandan . 2001 : Deegan et al.1. 12OO3I:" lYoneyama et al 
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Yoneyama et al.l . l2008l ) since the work of lYuse and Sand (|l993l ). In a 



typical experiment, a glass strip with a notch at its end is pulled at a constant velocity from 
an oven into a cold bath (see Fig. [T]). The control parameters in this experiment are mainly 
the width of the strip, the temperature gradient between the oven and the cold bath, and the 
pulling velocity. If the velocity is small enough, a crack does not propagate. Above a first 
critical velocity, the crack starts propagating following a straight centered path, and above a 
second critical velocity, the crack begins to oscillate with a well defined wavelength. However, 
the nature of the transiti on from straight to wavy path is not fully understood. For example. 



the experimental results (Ronsin, 19961 : Yuse and Sanol . 199?! ) are not sufficient to determine 
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whether the bifurcation is super-critical or sub-critical (also referred to as continuous or dis- 
continuous transition|3- At higher pulling velocities, irregular oscillations and branching are 
observed. Interestingly, these are very reproducible regimes, which makes the thermal crack an 
ideal model experiment as it allows to study the slow propagation of cracks under well controlled 
conditions. 



In an attempt to provide a theoretical explanation, iMardeii (|1994| ) successfully determined 
the propagation threshold by calculating Ki for a straight crack propagating at the center of the 
plate. He also proposed a descripti on of the oscillatory instability using the stability criterion 
derived by Cotterell and Rice ( 1980l ) for a crack in an infinite plate (the celebrated T-criterion), 
but that implied the improbable requirement that the fracture energy should depend strongly 
on the ve locity. Event u ally, this result was found to be incompatible with the experimental 



evidence ( Ronsin et al. , 19951 ). and serves as a classical example for the violation of the T- 



criterion. ISasa et al.l (| 19941 ) attempted to predict the threshold of the oscillatory instability and 
its wavelength by using the PLS under the infinite-plate approximation. Although they gave a 
reasonable qualitative description, their re sults deviated systematically fro m the experimental 
measurements. The next step was taken by Adda-Bedia and Pomeau ( 19951 ). who calculated Kn 
directly for a sinusoidal perturbation in a finite strip, in the configuration where the crack tip is 
at the center of the strip. However, in order to predict the oscillatory instability they resorted 
to an additional criterion, namely that the stability w i th res pect to this perturbation depends 
on the sign of Kn. More recently, iBouchbinder et al.l (120031) wrote an amplitude equat ion for 
a sinusoidal perturbation using the propagation criterion of Hodgdon and Sethna ( 19931 ). One 
should note that all these approaches rely on additional criteria when determining the instability 
threshold and cannot do so using only the classical principles of crac k propagation. In addition. 



theoretical results, except an attempt in (jBouchbinder et al.l . l2003l ^. are mostly limited to the 
study of the linear stability of a straight crack or to the study of transient regimes. This 
situation calls for a reexamination of the thermal crack problem in order to explain it using 
solid arguments. 

In this work, we begin with the study of quasistatic cracks in a quenche d glass plate using 
a pha se -field approach. P h ase fi eld models were originally introduced by ^ Caginalp and Fife 
( 19861 ): Collins and Levine ( 1986 ) to describe the propagation of solidi fication fronts. Ex- 
tensi ons of this work allowed quantitative modeling of dendritic growth (jKarma and Rappel . 
19981 ) and since then, phase field n iethods have been appl ied to many solidification prob- 



lems s uch as growth of binary a l lovs (lEtchebarria et al.l . 12004 ) and formation of polycristalline 
solids (jKobayashia and Warrenl . l2005l ). Furt hermore, they hav e been extended to ot her free 



boundary problems including viscou s fingering (iFolch et al.l.l200d) and crack propagation (iKarma et al. 
200ll : iKarma and Lobkovskvl . l2004l : iHenrv and Levinel . l2004l : iHakim and Karmal . l2005l . l2009l )" 



In the latter case, the phase field approach turns out to be very similar to the variationa l 
approach of fracture whic h uses the description of free discontinuity problems (Braidei. 19981 : 
Acerbi and Braides . 19991 ) using the F-convergence dAmbrosio and Tortorelli Il990l . Il992 ^. The 



numerical results of the phase-field model we derive allow to clarify the nature of the transition 
between straight and oscillatory cracks. 

We then pr ovide a theoretical analysis of the experiment of lYuse and San 3 (|l993l ). As i 



previous works (l Adda-Bedia and Pomeaul . Il995l : IBouchbinder et al.l . l2003l ) , the analysis is based 



m 



on the calculation of the SIFs. However, we calculate them to first order with respect to a straight 
centered crack for an arbitrary trajectory in a strip of finite width. We are then able to determine 
the instability threshold and the wavelength of the oscillations using only the PLS without 



^In (|Yuse and Sand . Il997l ). the bifurcation is claimed to be supercritical. Nonetheless, dat a points close to 
thres hold are not accurate e noush to rule out any other hypothesis (see Fig. 5 on page 372 in (|Yuse and Sand . 
Il997l )'). Results presented in l|Ronsiril Il99d ) suffer from the same lack of data points close to threshold. 
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Figure 1: The experimental setup. A glass plate of width 26 is moved slowly at a constant 
velocity v from an oven maintained at temperature Tq + AT to a cold bath of fixed temperature 
Tq. The crack tip is steady in the laboratory frame (moving in the plate frame) and its position 
with respect to the cold bath (bl) is controlled by the temperature field. 



introducing any additional criterion. Results agree very well with t he numerical simulations 
using the phase field model. And comparing them with the results of lAdda-Bedia and Pomeau 
(I1995I I. we find a small correction for the i nstab ility threshold, and a significant one for the 
wavelength. Also, unlike Bouchbinder et al. ( 20031 ). we argue that the calculation of the SIFs can 
be combined with the standard crack propagation criteria (that is with the Irwin criterion and 
the PLS) to obtain an evolution equation for the crack tip. Our main conclusion is that the PLS 
provides a good description of crack paths without any additional criteria. Recen tly, a similar 



appr oach for the thermal crack problem has been proposed and solved numerically (jPham et al 



20081^ which is in agreemen t with both our theoretical findings and with experimental results 
of Ronsin and Perrin ( 19981 ). 



The paper is organized as follows. We first provide a description of the thermal crack 
problem. We then describe the phase field approach for this problem, and compare its results 
with those of a theoretical analysis we perform. An important result of the theoretical analysis 
is an evolution equation for the crack path which allows quantitative predictions. We conclude 
by discussing the relevance of this work to path prediction of cracks in general. 



2 The thermal crack problem 

We first introduce the thermal crack problem and set the basic definitions and notations. Here, 
we consider the idealized problem of an infinitely long strip containing a semi-infinite crack. The 
coordinate system is chosen so that the axis of symmetry of the strip corresponds to y = and 
the tip of the crack lies at x = (see Fig. [T]). We assume that the temperature field is uniform 
across the thickness of the strip, and that the strip is under plane stress conditions, so that the 
problem is actually two-dimensional. We will also use dimensionless variables, measuring the 
lengths in units of the half-width of the strip b, the temperature in units of AT, the strains in 
units of aT^T/b, and the stresses in units of Eax^T, where E is the Young's modulus and 
aT the coefficient of thermal expansion. As will be shown, the behavior of the system is mainly 
governed by two dimensionless parameters : the Peclet number P = bv/D, where D is the 
thermal diffusion coefficient and the dimensionless material toughness Kic = Kic/{EaT^T^/b). 

Provided the velocity of the plate is not too low, the temperature field, which is assumed not 
to be affected by the presence of the crack, can be described by the advection-diffusion equation 
dxxT + PdxT = 0, with the boundary conditions: T{x) = in the cold bath and T{x) = 1 for 
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oo (jMardeij . Il994l ) . The solution of this equation is 



Tix) 



-P(x+l) 



@{x + l) 



(4) 



where x = —I corresponds to the location of the surface of the cold bath (see Fig. [T]) and Q{x) is 
the Heaviside function. Recall that the temperature field as given by Eq. ^ is an approximation 
of the experimental one. For a quantita tive comparison with experiments, the real temperature 
field should be used ( Pham et"all . boosi ). 

We consider the problem of a crack propagating through the heated strip within the frame- 
work of linear elastic fracture mechanics. Under plane stress conditions, the two-dimensional 
stress tensor a is related to the two-dimensional strain tensor e by 



1 



(J, 



1 



1 - v)eij + vekAj - (1 + u)T{x)5ij\ 



(5) 



where v is the Poisson ratio and the strain tensor e is related to the displacement field u by 



dui duj 
dxi dxj 



(6) 



The strip is free from external traction, thus within any approach one should always satisfy the 
boundary conditions 

ayy{x,±l) = axyix,±l) = . (7) 



3 The phase field model 



The classical theory of crack propagation, where the crack behavior is determined by the singu- 
larities of the stress field at the crack tip leads to difficult numerical issues when considering the 
movement and interaction of many cracks or to track crack motion in three dimensional systems 
where even the dynamics of crack f ronts is not we l l understood. In this context, the phase 
field approach to crack propag ation ( Aranson et al. . 200d : Karma et al. . 2001 : East gate et al 



2OO2I : iMarconi and .Taglal . Wm is very useful as it allows to study problems in arbitrary geome- 
tries and go beyond simple crack paths. It has succeeded in reproducing qualitative behavior 
of cracks such as branchi ng (jKarma and Lobkovskvl . 12004 ) and oscillations under biaxial strain 
( Henrv and Leving . 2004] ) . More recently, it has been shown that the Griffith criterion and the 
PLS are embedded in the phase field model of crack propagation ( Hakim and Karma . 20091 ). 



The general idea of phase field modeling is to introduce an additional field (the phase field) 
that describes the state of the system. The main advantage of this approach is that one does not 
need to treat the interface explicitly since it is defined implicitly as an isosurface of the phase 
field. This advantage becomes clear when studying the evolution of complex shapes since the 
algorithmic cost of using an additional field is much smaller than the cost of dealing explicitly 
with complex moving surfaces. In the case of fracture, this phase field 4> indicates whether the 
material is intact ((/> = 1) or broken (0 = 0). In the usual sharp interface representation, a 
crack surface is an infinitely thin boundary between a region where the material is intact and an 
empty region that can not sustain any stress. In contrast, the equations of the phase field model 
are such that (p varies continuously in space, and a crack surface is represented by a region of 
finite thickness. The thickness can be chosen freely and does not affect numerical results as long 
as it is much smaller than the characteristic length scale of the studied phenomenon. It can be 
shown that through an appropriate coupling between the phase field (p and the elastic fields, 
the model can have the desired properties : no stresses are transmitted across a crack (in the 
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limit where the system size is much larger than the width of the diffuse interface) and one can 
associate a finite surface energy with the interface, so that the fracture energy is well-defined. 

Here, we use this model to describe a quasistatic crack propagating due to a thermal gra- 
dient. We begin by presenting the model, including the adjustme nts needed to incorporate 
therm oelastic effects. The dimensionless elastic energy is written as (jAdda-Bedia and Pomeaul . 
I995I I 



E. 



el 



dxdy g{ 



1 



4(1 



1 



2{l + iy) 



(8) 



where the function g{(j)) = 0^(4— 3(/>) describes the coupling between the phase field and the elas- 
tic fields. This function is such that in regions where the material is broken ((p = 0), the contri- 
bution to the elastic energy is zero, while in regions where the material i s intact, the contrib ution 



to the elastic energy recovers the one prescribed by linear elasticity. In (jKarma et al.l . 120011 ). it is 
shown that the particular choice of g does not affect the results as long as g{0) = g'{0) = g'{l) = 
and Iim0=o5(</') ~ 4'"^ with a > 2. The additional term T{x) corresponds to the thermal ex- 
pansion. This effect is assumed to be the same in the intact and in the partially broken parts 
of the material. 

Since we are considering a quasi-static crack growth, in which the elastic body is at me- 
chanical equilibrium at all times, the elastic energy should be an extremum with respect to the 
displacement field u{x,y). That is 

^ = . 9 

dUi 



We now need to specify the evolution of the ph ase field. This is done by associating with the 
phase field a dimensionless free energy given by (|Henrv and Levind . |2004| ) 



E, 



where £s is defined by 



dxdy 



(10) 



£<h 



-h) i^kk - 2T{x)f + 2(1^ (eij 



1 



en 



-^kk 



if {eu - 2T{x)) > 
if {eu - 2T{x)) < 
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where V{cj)) = 1.5/i(l — 0^)(/>^ is a double well potential, with an energy barrier of height oc h, that 
ensures that the preferred states of the homogeneous system are either = 1 (intact) or = 
(completely broken). £^ coincides with the standard elastic energy density when the material 
is under tension and includes only the contribution of shear when the mat erial is under com- 
press ion. This choice avoids the propagation of a crack under compression (iHenry and Levind . 
20041 ). Taking into account the coupling through the function g{<j3), when f"^ is larger than a 



threshold value £c, the broken phase is favored while when £^ < £c, the intact phase is favored. 
The evolution equation for (p is relaxational with a kinetic bias that ensures that the material 
will not be driven from a broken state to an intact state : 



mm 



5E^, 



,0 



:i2) 



where r is a constant. The irreversibility of Eq. (jl2p means that a crack cannot heal (which 
would happen in our set up at the crack tail) and ensures that the total elastic energy stored in 
the material cannot increase when 6 varies. 



6 



The elastic equation ([9]) is solved on a rectangular domain [— x [—1, 1] w ith L ^ 1 using 



the Gauss-Seidel over-relaxation method at each time step (jPress et al.l . 120071 ) . In addition to 
the boundary conditions d?]), one has to set additional boundary conditions at x = ibL/2. Here 
we have chosen 

eiji-L/2,y) = 0, (13) 
a^y{L/2,y) = a^^{L/2,y) = . (14) 

It is clear that the final solution is insensitive to these boundary conditi ons as long as L ^ 1 /P. 



The phase field equation (jl2p is solved using an Euler forward scheme (jPress et al.1 . 120071 ) . The 
values of the grid spacing dx and of the diffusion coefficient are chosen so that the diffuse 
interface occupies few grid points and is much smaller than the characteristic wavelength of 
the oscillating crack, while the characteristic time r is chosen so that the phase-field relaxes 
fast enough to its equilibrium (or more precisely, multiplying r by a factor 2 does not affect 
quantitatively the results.). Numerical checks were performed to ensure that neither changes in 
the convergence criterion of the Gauss-Seidel iterative method nor changes in the grid resolution 
affect the results. The size of the domain [L] x [2] was typically 1500 x 300 grid points, and 
numerical checks showed that doubling L did not a,ffect the resu l ts. Fi nally, the fracture energy 
r = K^^/ (2/i) was computed using the expression ( Karma et al.l . I2OO1I ) 



r = d<p^2D4£,-i£,g{<l))-V{cP)). (15) 
leading to the dimensionless material toughness Kjc = \/2fiT /{EaT^T\/T)). 



Before turning to the description of numerical results, we find it necessary to summarize the 
different parameters used here and to give some rationale for their choice. In addition to the 
physical parameters associated to the material and to the temperature field, the pertinent phase 
field parameters are r, D^, h and Sc- In dimensionless units the values of r = 1, = 4.4 x 10^'^, 
h = 1 and £c = ^ have been chosen. Notice that the three latter parameters determine both 
the fracture energy and the interface width of the phase field. Since we are considering the 
dimensionless material toughness Kjc, the only effect we are concerned with, is the dependance 
of the width of the diffuse interface (the broken surface) on Z)^. More precisely, in the spirit 
of the phase field approach, the value of D^, for a given value of Kjc should not affect the 
behaviour of the crack. Within the present simulations, we have checked that using a value of 

corresponding to an interface which is twice thinner does not affect quantitatively the final 
results. 

We now summarize the results of the simulations. By tuning the control parameters P and 
Kic, we could observe the different regimes observed in experiments : no crack propagation, 
propagation of a straight crack at the center of the strip or a crack following a wavy path (see 
Fig. [2]). However, the main purpose of the present simulations was to determine the nature of the 
bifurcation from straight to oscillatory cracks, since no direct experimental results are available. 
Namely in experiments, the amplitude of oscillations varies almost linearly with the control 
parameter, but since experimental results do not allow to explore the behavior of the crack 
close to the threshold it was impossible to characterize the crack path in the close vicinity of the 
transition point. Furthermore, it should be em phasized that apa r t from attempts describin g this 



transition by a Landau-Ginzburg expansion (|Bahr et al.l . Il995l : iBouchbinder et al.l . l2003l l. the 



available theoretical work has been mostly limited to linear stability analysis. Using the phase 
field model and without resorting to any assumption, we were able to compute the steady-state 
paths and to determine the nature of the bifurcation. 

First, in Fig. [3] we present some typical crack paths close to the threshold for the sake of 
comparison with the analytical results presented below and also in order to clarify the nature 
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Figure 2: Typical crack paths obtained using the phase field model. The white region corre- 
sponds to the broken surface and the color code corresponds to the density of elastic energy in 
the material (red: high elastic energy density, blue: zero elastic energy), (a) A straight centered 
propagating crack, (b) A wavy crack above the instability threshold. This picture shows that 
the crack width is much smaller than the wavelength of the observed oscillations. 
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Figure 3: (a): Crack tip trajectory below the threshold of linear instability for a crack which 
was initially off center. P = 7.9 and 1/Kic = 7.85. (b): Crack tip trajectory for an initially 
slightly off center crack above the threshold of linear instability. P = 7.9 and 1/Kic = 8.43. 



of the bifurcation. The first path was obtained with an initial straight crack which was off 
center and below the threshold. One can see damped oscillations as the crack path returns to 
its equilibrium state: a straight crack. The second path was obtained with an initial crack that 
was slightly off center (shifted by 2 grid points away from the center) and slightly above the 
threshold. One can see oscillations of the path that are amplified as the crack advances. These 
oscillations eventually saturate and the crack path becomes oscillating with a finite constant 
amplitude as observed in experiments. Using the model we were able to compute the amplitude 
of the oscillations when a steady state was reached. The simulations showed that the steady 
state is independant of initial conditions, such as the initial position of the crack, and is a single 
valued function of the control parameters. Typical results are shown in Fig. [4] where the square 
of the amplitude is plotted as a function of 1/Kic for a given value of P. One can see that 
above a threshold value of 1/Kic the amplitude of oscillations is no longer zero and scales as the 
square root of the deviation from the threshold. The same behavior holds when using the value 
of P as a control parameter. The amplitude curves together with the damped oscillations below 
the threshold and the amplified oscillations above the threshold indicate that the bifurcation is 
supercritical or continuous. 

For values of the control parameters well above the thr eshold, we wer e not able to re- 
cover complex paths similar to those observed in experiments (lYuse and San I I1993I I. This may 



be due to the fact that well above the threshold the quasistatic approximation is no longer 
valid and dynamic effects through the temperature field as given by Eq. [4] become relevant. 
It is nonetheless interesting to note that the observed crescent-like path in Fig. [5] exhibits a 
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Figure 4: (a) Square of the amplitude A of the oscillations for P = 7.89 as function of 1/Kic, 
which is related to the amplitude of the temperature gradient, (b): same as (a) with P = 3.95. 
The linear behavior of A"^ in the neighborhood of the transition shows that the transition from 
straight to oscillatory propagation is supercritical (continuous). 
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Figure 5: Crack tip trajectory for a crack deep in the nonlinear regime, that is far above the 
threshold of the instability. The values of the control parame ters are P = 7.9 and i/K]c 



11.6. Note the strik i ng resemblance to crack p aths appearing in (jGhatak and Mahadevanl . 12003 



Audolv et all l200,4 ISendova and Willisl . [2003) 
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Figure 6: Main theoretical and numerical results, (a) The phase diagram of the different states 
of the system as a function of the two control parameters P and K\c. (b) The wavelength A 
as a function of P at the transition from straight to wavy crack propagation. In both figures, 
the symbols + and x correspond to the numerical results obtained from the phase field model, 
the solid lines correspond to the anal ytical calculation presented in Se c. U] and the dashed lines 
correspond to the results reported in ( Adda-Bedia and Pomeau . 19951 ). 



strikin g similitude with paths obtained when cutting a th in elastic sheet with a moving thick 
object ( Ghatak and Mahadevan . 20031 : Audolv et al. . 20051 ) and with crack paths obtained dur- 
ing drying of silicate sol-gel films fabricated using spin-coating techniques ( Sendova and Willis! . 

Eos). 

The phase diagram in Fig. [6l^a) shows the threshold for the propagation of a centered straight 
crack and for the transition to a wavy crack propagation in the l/P-\/Kic phase space. It 
can be seen that the theoretical predictions presented in Sec. H] are in quantitative agreement 
with numerical results of the phase field model without any adjustable parameters. The same 
agreement is reached when considering the dimensionless wavelength at the threshold of the 
transition from straight to wavy crack path (see FiglGj^b)). For this latter quantity, one should 
also note the good agree nient with experimen tal values for small values of 1/P (linear behavior 
with A ~ 0.28 + 2.0/P) dRonsin et al.l . Il995l '). Moreover, the phase field computation and our 
theoretical approach lead to the same deviation from the linear beh avior for 1/P > 0.2, where 
the d iscrepancy between the present results and those reported in ( Adda-Bedia and Pomeaul . 
19951 ) become significant. This feature has also been observed in experiments (jRonsinl . Il996l ). 



4 Theoretical analysis 



We now turn to a theoretical study of the tran sition from straight to oscillatory crack propa- 
gation. As in previous t heoretical treatments (IMardeii . Il994l : lAdda-Bedia and Pomeaul . Il995 



Bouchbinder et al.l . l2003l ^. the analysis is based on the calculation of the stress intensity factors 



(SIFs). However, here, they are computed to first order for an arbitrary trajectory without 
using the infinite plate approximation and without the restriction that the crack tip lies at the 
center of the sample. Within the framework of linear elastic fracture mechanics, the propagation 
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of a crack is governed by the singular behavior of the stress field in the vicinity of its tip, as 
defined in Eq. ([T]). The propagation laws we use are the Irwin criterion defined in Eq. ([2]) and 
the principle of local symmetry (PLS) defined in Eq. ([3|). To calculate the SIFs, one must solve 
the equilibrium equations 

, vVii = -v2r(x) , (16) 



with the boundary conditions ([7]) at the sides of the plate and the additional boundary condition 
at the crack surface 

aijUj = , (17) 

where ft is the normal to the crack faces. 

To investigate the transition from a straight to an oscillating regime, we consider a small 
deviation from a straight centered crack. 



y{x) = Ah{x) 



(18) 



where h{x) is defined for < (see Fig. [T|) and \A\ <^ 1. We expand the displacement and 
stress fields to first order as 



a. 



(0) 

i 

M 



+ Auf'> + O {A^) , 



(19) 
(20) 



The SIF Ki (resp. Kn) is an even (resp. odd) function of A. Therefore their expansions have 
the form 



Kii 



k(0) + o{A^) 
Ak[1^ + 0{A^ 



(21) 
(22) 



Note that in particular, Kn = for a straight centered crack. The different terms of the above 
expansion can be determined by solving the equilibrium equations using repeatedly the Wiener- 
Hopf method. The details of these calculations are presented in Appendix A. Here we simply 



summarize the calculation of , which is novel. can be written as 



K^\h{x)] = Kh{0) + Kl^'iHx) - h{0)] , 



(i)r 



(23) 



which amounts to decomposing the general problem into that of a straight off- center crack and 

that o f an oscillating crack with a centered tip. The latter problem was solved bv lAdda-Bedia and Pomeau 
(|l995l l. yielding 



KS\h{x)-hm 



(0) 



-/i'(0) + 



dx 



d_ 

dx V 



a^J>Jix,0)ihix)-hm 



p' 



(24) 



where p^{x) is a weig ht function that does not d e pend on the physical parameters. We have gen- 
eralized the results of Adda-Bedia and Pomeau (|l995l ) with the calculation of the contribution 
Kh{0) to obtain 



4'>[/i(x)] 



a 



(0) 



K 



(0) 



d 



HO) + -^h'{0) + / dx— [4°j(x,0)/i(x) 



(25) 



where c ~ 1.272 is a numerical constant (see Appendix A for details). 

Having calculated the SIFs to first order, we can now examine the particular form taken by 
the crack propagation criteria for a small perturbation. The Irwin criterion reads K^^^ = Kic. 
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t igure as a function of I, for P = 5. Recall that / is the distance of the crack tip from 

the boundary of the cold bath. The threshold of incipient growth of a centered straight crack is 



given by Ki^ = K^^ (/, 



where Ic is the location of the maximum of K-^ 



(0) 



{I). 



The transition to a propagating straight crack is thus governed by the existence of I such that 



Kf\l) = Kic. This occurs when Kic = (^c); where Ic is the location of the maximum 
of K^^\l) (see Fig. [7]). This condition determines the transition curve in Fig. [2] from no crack 
propagation to a propagation of a centered straight crack. Also, above the propagation threshold 
Eq. ()2ip shows that to first order in A, the distance between the cold front and the crack tip, 
which is determined by Kj [h{x)], does not depend on the crack path y{x) = Ah{x). This 
simplifies the problem, as we can parameterize the "time" evolution of the crack path by the 
position of its tip in a coordinate system attached to the plate. If h{x) now describes the crack 
path in such a coordinate system, the value of K^^^ when the tip is at the point x is simply 
given by 



'(0) 



K^\x) 



(0) 



Kx) + —k-h'i^) + I du 



d_ r 

du 



a. 



(0)/ 



{u,0)h{x + u) p^{-u) . (26) 



Therefore, the PLS which for a small perturbation reads K^j^ (x) = 0, yields 



cK, 



(0) 



d_ 

dl 



(0) 



h{x) + -^h'ix) + 



Q 

du 

-oo 



du 



a^^]{u,0)h{x + u) p^{-u). (27) 



Eq. (j27p is a linear equation for a crack that deviates slightly from the straight centered path. 

To determine the stability of the straight centered crack propagation, we must examine the 
evolution of a small perturbation, which amounts to solving the following problem : given a crack 
path h{x) defined for x < 0, which does not necessarily satisfy {x) = for x < 0, find a 
continuation of h such that K^\x) = for x > 0. Now, since cj^x{x,^) decreases exponentially 
as X — > — oo, one can expect the long-term behavior of the crack to be the same as that of a 
crack path satisfying {x) = for all x, and this is indeed observed in the simulations (see 
Fig. [3Kb) in the previous section, and the discussion below). Therefore, examining crack paths 
that satisfy K^^^ = from the beginning is sufficient to determine the stability threshold, and 
we will use such paths for that purpose. 
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We look for solutions of the form 



h{x) = Re [6 eKp{iqx)] 



(28) 



where Re stands for the real part, 6 and q are in general complex numbers. Then using Eq. (I26p . 
K^^\x) takes the form 

K^l^ (x) = Re [5 Kuiq) exp(iQx)] , (29) 

where 



Kuiq) 



A solution K^^\x) = is obtained if Kii[q) = 0, which is the dispersion relation of the problem. 
Since q and Kii{q) are in general complex numbers, we can expect the equation Ku{q) = 
to admit a discrete set of solutions for a given set of the dimensionless parameters P and Ki^. 
By examining the behavior of Kii{q), we find that the stability is determined by solutions 
corresponding to a pair of conjugate oscillating modes, which are shown on Fig. [SKa) for a given 
value of P. When lm(q) > 0, the perturbation (j28p is damped which implies that the straight 
and centered crack propagation is stable. When lm(q) < 0, the perturbation is amplified, leading 
to a wavy crack propagation regime. Thus the instability threshold corresponds to Im(g) = 0, 
and is denoted hy q = qc- Using this result, another way to determine the instability threshold is 
illustrated in Fig. [81(b). It consists in supposing from the beginning that q = u is a real numbeiEl 
and in looking for the specific value lo = qc for which Im[Jrn(w)] = and Re[i^n('^)] = 
simultaneously. 

By repeating the above computation for different values of P and Kic, we obtained in Fig.[6l^a) 
the phase diagram and in Fig. M^h) the wavelength of the oscillations at the threshold, which 
is given by A = 2Tr /qc- These figures show that a good quantitative agreement is reached with 
the numerical results of the phase field model for both the wavy instability threshold and the 
wavelength at the transition. 

Our results for the instability threshold differ from those of Adda-Bedia and Pomeau ( 19951 ). 
in which the threshold was determined by the existence of a real solution qc such that K^j"^ [sin lux] = 

^ [sin ujx] = for x = and the stability criterion was based on the sign of -fCjj [sinwx]. In 

general, such a solution has no reason to satisfy K^'' (x) = for x 7^ . Indeed, the wavelengths 
predic ted by these two criteria (i.e., the present one and the one of lAdda-Bedia and Pomeau 
differ significantly, especially for small values of P (see Fig. [6)^b)). The thresholds pre- 
dicted also differ, but less notably (see Fig. El^a)). As can be seen in Fig. [8l[b), the deviations 
between the estimates of Ki^ and A obtained by the two methods are related to the vertical and 
horizontal distances between the minimum of Im(i^n) and its intersection with Re(-ftrn)) and 
the vertical distance scales as the square of the horizontal one. 



(0) 



d_ 

dl 



(0) 



+ iq- 



(0) 



+ 



Q 

du— 

^nr, OU 



cT£)(n,0)e 



iqu 



(30) 



5 Conclusions 

In this paper we address the problem of crack path prediction for quasi-static crack propagation. 
We used the framework of the thermal crack experiment since it is a well controlled experiment, 
and the best known to us. The reason is that stresses are induced internally by an imposed 

^we use LO instead of q in order to emphasize the fact that it is a real number 



13 



0.15 




I 1 i i 1 _o.i5 I , , , 1 

8 9 10 11 12 0.5 1 1.5 2 

l/iTic t^/27r 

Figure 8: (a) Solutions of Kii(q) = versus l/i^ic for P = 5. The solution branch disappears 
when lm{q) exceeds a certain value, as the integral appearing in the calculation of Kii{q) ceases 
to converge, (b) The real and imaginary part of KuIlu), with uj real, at the instability threshold 
for P = 5 {Kic ~ 0.103). 



thermal gradient that can be easily tuned by external parameters leading to a steady state 
propagation. 

Our main result is that by using only the Principle of Local Symmetry and the Griffith 
criterion one can reproduce the phase diagram of the problem, namely predict the existence 
of the two instabilities (no-crack to straight crack and straight to oscillatory propagation). 
Furthermore, we show that these criteria are sufficient to determine the instability threshold 
and the wavelength of the oscillations without introducing additional hypotheses. Interestingly, 
our theoretical results agree very well with the numerical simulations of the phase field model. 
This strengthens the theoretical approach based on the classical LEFM theory combined with the 
two propagation criteria, namely PLS and Griffith criterion (but mostly PLS). This work joins 
a recent effort to exp lain a variety of instabilities of cracks, namely the branching in s tability 
(|Katzav et al.l . l2007bl l and the out-of-the-plane roughening instability (jKatzav et al.l . l2007al ) 
using these general criteria (i.e., LEFM + Griffith criterion + PLS), 

Another achievement of this work is the fact that the phase field model allows to clarify the 
nature of the transition between straight and oscillatory cracks, which is convincingly shown 
to be supercritical or continuous. This is achieved thanks to the high resolution of the phase 
field approach near the transition point and the relative ease at which in-silico experiments 
can be performed under different conditions. This also suggests that the phase field model is 
more adapted than the classical methods for crack path tracking, especially when considering 
complex moving shapes, in view of its relatively low computational cost. It is therefore a good 
can didate to describe complic ated crack pa tterns such as s pirals or crescents as t he ones observed 



in (ISendova and Willid. 120031). br anching (jKatzav et al 



2007bl : iHenrv I . I2OO8I ) and interacting 



cracks ([Marconi and Jaglal . l2005l ). A real challenge for this approach would be to reproduce 
propagation of cracks in 3D, which can help to infer their laws of motion from such well controlled 
numerical experiments. 
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A The straight off-center crack 



In the following, we present briefly the method of resolution for a stra ight off-center crack in a 
therm al gradient. The method of resolution follows the same steps as in (jAdda-Bedia and Pomeaul . 
19951 ). Consider a straight crack which is slightly above/below the center of the strip, namely 
its location is given by y{x) = Ah{0) = 5. Due to the absence of symmetry with respect to the 
center of the strip y = 0, let us split the displacement, strain and stress fields in the strip as 
follows 

f{x,y) = ^{f+{x,y)&{y - 5) + r{x,y)ei5 - y)) , (31) 

where Q is the heaviside function and / stands for any component of the elastic fields. Let us 
also define 

[f]ix) = ^{f+{x,6)-r{x,6)) , (32) 

Note that all the elastic fields are continuous in the unbroken regions. However, along the crack 
surfaces the displacements might be discontinuous while ayy and axy are continuous there. Thus, 
for this specific problem, the boundary conditions (j7ll7p can be rewritten as 



a^y{x,±l) 
(x) 



aty{x,±l)=0, 



Cf. 



yyi 



(7: 



xy 



a, 



[Ux] (x) 

^'x,6) 



yy 



a, 



xy 



] (^) = , 
{x) = 
{x,5) = 



for X > , 
for X < 



(33) 
(34) 
(35) 
(36) 



Now, we Fourier transform along the x direction all the elastic fields and the temperature 
field and plug them into the equilibrium equations (fT6]) and into the boundary conditions (j33l34p . 
After some algebraic manipulations, this yields relations of the type 



[Uy] (k) 



Cxx{k) Cxy{k) Dx{k) 

Cyx{k) Cyy{k) Dyijz) 



axy{k,6) 
ayy{k,5) 

f{k) 



(37) 



Performing the first-order expansions in 5, as given by Eqs. (jl9p -()20p: 

Ui{k,5) = uf\k,{)) + 5uf\k,{)) , 
^^Jik,S) = alfik,0)+6aff{k,0), 

and inverting (|37p gives the final result 



(38) 
(39) 



where 



<7g)(A:,0) 
a£)(A:,0) 
4S(fc,0) 



u. 



(0) 



{k,{))+D{k)f{k) 



H{k)a^^^{k,0) + S{k)T{k) , 



F{k) = k 



-P{k) 
sinh^ k — k"^ 



H{k) 
P{k) = k 



sinh2A: + 2k 
k'^ + sinh^ k 

sinh^ k — k"^ ' 
sinh^ k — k"^ 

sinh2/c - 2k 



Dik) 
S{k)-. 



(40) 
(41) 
(42) 



2 (1 — cosh A;)(sinh k — k) 

sinh 2k + 2k 
sinh k — k 



sinh k + k 



(43) 
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The first step of resolution wliicli involves only the zeroth-order expansion, i.e. the case of a 
straight centered crack, is the calculation of K^^^ . Using the boundary c onditions (I35I36D to lead- 
ing order and applying the Wiener-Hopf method to equation ()40p yields (jAdda-Bedia and Pomeaul . 
\l995t ) 



(0) 



+00 



D{k)T{k)F+{k)- , 



where 



F-{k) 



F+{k) ' 

and F^{k) {F^{k)) is analytic in the upper (lower) half-plane, respectively. 



(44) 



(45) 



To complete the calculation of the SIFs, we return to the case of a straight off-center crack 
and determine the coefficient k appearing in Eq. ([23l) . Using the boundary conditions (|35l36p 
to first order in 5 and applying again the Wiener-Hopf method to Eq. (I42p yields 



u. 



(1) 



{k) 



ik 



4°' 



k 



ik) 



1 



p-{k) y_ 



ik'P+{k')af}~{k',{)) dk' 



+ 



a 



k'-k + ie 2iTT P-{k) 



F-{k) 
1 



D{k')f{k')F+{k')dk' 



-co 

+00 



k + ie 27r 

{k',0) dk' 



k' -k 



+ 



a 



2m P-{k) ' 



(46) 



where P{h) = P {k)/P^(k) with P^{k) (P (k)) analytic in the upper (lower) half-plane, 
and e is a vanishingly small positive number, notice that the solution of u^^ (k) obtained 

from the zeroth order calculations has been used, and that <t^x{x,0) can be obtained from 
Eq. (|4ip. For both F{k) and P{k) the factorization i nto F^(k) and P^(k) can be achieved 



using a Fade decomposition, as previously explained in ([Bouchbinder et al.l . l2003l : iKatzav et al 



20073). The unknown constant a is introduced because the decomposition involved in the 

The value of a can be determined by noting that the 
00. Canceling out the k~^^'^ term on the r.h.s. yields 



Wiener-Hopf method is not unique. 
(k) behaves as 0(fe~^/^) as k - 



a 



dk 

D{k)T{k)F+{k)- 



(47) 



and thus one has 



{k) 



1 



1 



+- 



F-{k) 
1 



D{k')f{k')F^ 



+~ -ik'D{k')f{k')F+{k') dk' 



dk' 
'2^ 



F-ik) 

Finally by looking at asymptotic behavior of 



k' — k + ie 2m 
(k) as A; — > 00 one has 



(48) 



cK- 



(0) 



d_ 

dx 



a^^]{x,0)p+{-x)dx 



(49) 



where p'^{x) is the inverse Fourier transform of P'^{k) and c is a numerical constant given by 



lim 

fc— >oo 



(ik) 2 

"Tf" 



1 



1 



F-{k) P-{k) 



1.272 



(50) 



Using this value of n in Eq. (p3]) yields Eq. ([25]). 



16 



References 

Acerbi E, Braides A (1999) Approximation of free-discontinuity problems by elliptic functional 
via F-convergence. Asymptot. Anal. 21:317-329 

Adda-Bedia M, Pomeau Y (1995) Crack instabilties of a heated glass strip. Phys Rev E 52:4105- 
4113 

Adda-Bedia M, Arias R, Ben Amar M, Lund F (1999) Generalized Griffith criterion for dynamic 
fracture and the stability of crack motion at high velocities. Phys Rev E 60:2366-2376 

Ambrosio L, Tortorelli VM (1990) Approximation of functionals depending on jumps by elliptic 
functional via F-convergence. Comm. Pure Appl. Math. 43:999-1036 

Ambrosio L, Tortorelli VM (1992) On the approximation of free discontinuity problems Boll. 
Un. Mat. Ital. B 6:105-123 

Amestoy M, Leblond JB (1992) Crack paths in plane situations-IL Detailed form of the expan- 
sion of the stress instensity factors. Int J Solids Struct 29:465-501 

Aranson IS, Kalatsky VA, Vonokur VM (2000) Continuum field description of crack propagation 
Phys Rev Lett 85:118-121 

Audoly B, Roman B, Reis R (2005) Cracks in brittle elastic plates: when geometry rules fracture 
paths. Phys Rev Lett 95:025502 

Bahr HA, Gerbatsch A, Bahr U, Weiss HJ (1995) Oscillatory instability in thermal cracking: a 
first-order phase-transition phenomenon. Phys Rev E 52:240-243 

Barenblatt G. and Cherepanov G (1961) On brittle cracks under longitudinal shear. PMM 
25:1110-1119 

Bilby BA, Cardew GE (1975) The crack with a kinked tip. Int J Frac 11:708-712 

Bouchbinder E, Hentschel HGE, Procaccia I (2003) Dynamical instabilities of quasistatic crack 
propagation under thermal stress. Phys Rev E 68:036601 

Braides A (1998) Approximation of free-discontinuity problems. Springer 

Broberg KB (1999) Cracks and Fracture. Academic Press, London 

Caginalp G, Fife P (1986) Phase Field Methods for Interfacial Boundaries Phys Rev B 33:7792- 
7794 

Collins JB, Levine H (1986) Diffuse interface model of diffusion-limited cristal growth Phys Rev 
B 31:6119-6122 

Cotterell B, Rice JR (1980) Slightly curved or kinked cracks. Int J Fracture 16:155-169 

Deegan RD, Chheda S, Patel L, Marder M, Swinncy HL, Kim J, dc Lozannc A (2003) Wavy 
and rough cracks in silicon. Phys Rev E 67:066209 

Eastgate LO, Sethna JP, Rauscher M, Cretegny T, Chen CS, Myers CR (2002) Fracture in 
mode I using a conserved phase-field model. Phys Rev E 65:036117 

Erdogan G, Sih GC (1963) On the crack extension in plates under plane loading and transverse 
shear. Journal of Basic Engineering 85:519-527 



17 



Etchebarria B, Folch R, Karma A, Plapp P (2004) Quantitative phase-field model of alloy 
solidification Phys Rev E 70:061604 

Fineberg J, Marder M (1999) Instability in dynamic fracture. Phys Rep 313:2-108 

Folch R, Casademunt J, Hernandez-Machado A (2000) Viscous fingering in liquid crystals: 
Anisotropy and morphological transitions Phys Rev E 61:6632-6638 

Freund LB (1990) Dynamic Fracture Mechanics. Cambridge University Press, Cambridge 

Ghatak A, Mahadevan L (2003) Crack street: the cycloidal wake of a cylinder ripping through 
a thin solid sheet. Phys Rev Lett 91:215507 

Gol'dstein RV, Salganik RL (1974) Brittle fracture of solids with arbitrary cracks. Int J Fracture 
10:507-523 

Griffith AA (1920) The Phenomenon of Rupture and Flow in Solid. Phil Trans R Soc Lond Ser 
A Int J Fracture 221:163-198 

Hakim V, Karma A (2005) Crack path prediction in anisotropic brittle materials. Phys Rev 
Lett 95:235501 

Hakim V, Karma A (2009) Laws of crack motion and phase-field models of fracture. J Mech 
Phys Solids 57:342-368 

Henry H, Levine H (2004) Dynamic instabilities of fracture under biaxial strain using a phase 
field model. Phys Rev Lett 93:105504 

Henry H (2008) Study of the branching instability using a phase field model of inplane crack 
propagation. EPL 83:16004 

Hodgdon JA, Sethna JP (1993) Derivation of a general three-dimensional crack-propagation 
law-A generalization of the principle of local symmetry. Phys Rev B 47:4831-4840 

Irwin GR (1957) Analysis of stresses and strains near the end of a crack traversing a plate. J 
Appl Mech 24:361-364 

Karma A, Rappel WJ (1998) Quantitative phase-field modeling of dendritic growth in two and 
three dimensions Phys Rev E 57:4323-4349 

Katzav E, Adda-Bedia M, Derrida B (2007a) Fracture surfaces of heterogeneous materials: a 
2D solvable model. EPL 78:46006 

Katzav E, Adda-Bedia M, Arias R (2007b) Theory of dynamic crack branching in brittle mate- 
rials. Int J Fracture 143:245-271 

Karma A, Kessler DA, Levine H (2001) Phase-field model of mode HI dynamic fracture. Phys 
Rev Lett 87:045501 

Karma A, Lobkovsky AE (2004) Unsteady crack motion and branching in a phase field model 
of brittle fracture. Phys Rev Lett 92:245510 

Kobayashia R, Warren JA (2005) Modeling the formation and dynamics of polycrystals in 3D 
Physica A 356:127-132 

Leblond JB (1989) Crack paths in plane situations- I. General form of the expansion of the 
stress instensity factors. Int J Solids Struct 25:1311-1325 



18 



Leblond JB (2003) Mecanique de la rupture fragile et ductile. Hermes Science Publications, 
Paris. 

Marconi VI, Jagla EA (2005) Diffuse interface approach to brittle fracture. Phys Rev E 71:036110 

Marder M (1994) Instability of a crack in a heated strip. Phys Rev E 49:51-54 

Pham VB, Bahr HA, Bahr U, Balke H, Weiss HJ (2008) Global bifurcation criterion for oscil- 
latory crack path instability. Phys Rev E 77:066114 

Press WH, Teukolsky SA, Vetterling WT, Flannery BP (2007) Numerical Recipes. Cambridge 
University Press, Cambridge 

Ronsin O, Heslot F, Perrin B (1995) Experimental study of quasistatic brittle crack propagation. 
Phys Rev Lett 75:2352-2355 

Ronsin O (1996) Etude experimentale de la propagation de fractures dirigees en milieu fragile. 
PhD thesis, Universite Paris VI 

Ronsin O, Perrin B (1998) Dynamics of quasistatic directional crack growth. Phys Rev E 
58:7878-7886 

Sakaue K, Yoneyama S, Kikuta H, Takashi M (2008) Evaluating crack tip stress field in a thin 
glass plate under thermal load. Engineering Fracture Mechanics 75:1015-1026 

Sasa S, Sekimoto K, Nakanishi H (1994) Oscillatory instability of crack propagations in quasi- 
static fracture. Phys Rev E 50:1733-1736 

Sendova M, Willis K (2003) Spiral and curved periodic crack patterns in sol-gel films. Appl 
Phys A 76:957-959 

Yang B, Ravi-Chandar K (2001) Crack path instabilities in a quenched glass plate. J Mech 
Phys SoHds 49:91-130 

Yoneyama S, Sakaue K, Kikuta H, Takashi M (2006) Instantaneous phase-stepping photoelastic- 
ity for the study of crack growth behaviour in a quenched thin glass plate. Meas Sci Technol 
17:3309-3316 

Yoneyama S, Sakaue K, Kikuta H, Takashi M (2008) Observation of Stress Field Around an 
Oscillating Crack Tip in a Quenched Thin Glass Plate. Experimental Mechanics 48:367-374 

Yuse A, Sano M (1993) Transition between crack patterns in quenched glass plates. Nature 
362:329-331 

Yuse A, Sano M (1997) Instabilities of quasi-static crack patterns in quenched glass plates. 
Physica D 108:365-378 



19 



